A Model for Granular Texture with Steric Exclusion 
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We propose a new method to characterize the geometrical texture of a granular packing at the 
particle scale including the steric hindrance effect. This method is based on the assumption of a 
maximum disorder (entropy) compatible both with strain-induced anisotropy of the contact network 
and steric exclusions. We show that the predicted statistics for the local configurations is in a fairly 
agreement with our numerical data. 



In modeling granular media from a microscopic point 
of view, one is naturally led to consider a fundamental 
question: how do the properties of the spatial organi- 
zation of particles control the behaviour at the macro- 
scopic scale? A major difficulty is that the effects arising 
from the geometrical frustration of particles prevail due 
to hard core exclusions, and thus the standard methods 
developped for multibody systems fail as such to provide 
a satisfactory description of many subtle structural fea- 
tures of a granular medium. 

The most basic effect of mutual exclusions is that a 
particle may be simultaneously touched by only a few 
neighbouring particles. The coordination number (the 
number of touching neighbours) cannot exceed 6 in a 2D 
assembly of particles of nearly the same size and 12 in a 
3D packing. These local environments fluctuate strongly 
in space both in terms of the coordination number and 
the angular positions of the neighbours. 

Another basic observation is that the relative angular 
positions of the particles i.e. the directions of the contact 
normals are not isotropically distributed fl0|. 

The adequate description of local environments includ- 
ing these geometrical features, i.e. a combination of 
excluded- volume constraints, disorder and anisotropy, is 
a crucial step in tracing back the failure and flow proper- 
ties of granular materials to the particle scale. Consider, 
for instance, a granular packing in statistical equilibrium. 
The force balance on a particle involves the angular po- 
sitions of the few neighbours belonging to the local envi- 
ronment. The fact that two neighbours may not occupy 
two angular positions separated by an angle smaller than 
a finite angle 59 (approximately -| for particles of nearly 
the same size, see figure 0) reduces drastically the acces- 
sible equilibrium states of the particle and hence that of 
the whole packing. This effect is even more critical at 
incipient failure and for larger anisotropies of the con- 
tact directions. While the texture anisotropy has been 
considered by several authors in the past jij, the key 
role of steric exclusions has been almost completely dis- 
regarded in micro-macro models of the quasi-static rhe- 
ology of granular materials. 



In this paper, we propose a model that allows to con- 
struct local environments incorporating both the steric 
exclusions and generic disorder of granular media from 
the global geometrical texture described in terms of the 
directions of interparticular contacts or the correspond- 
ing anisotropy. The basic assumption of our model is a 
maximum directional disorder consistent both with the 
texture anisotropy and with steric constraints. We ap- 
ply our approach to a two-dimensional medium and solve 
the resulting equation. The predicted statistics of local 
environments will then be compared with raw data from 
numerical simulations. 

Let us consider the simple model of a 2D granular as- 
sembly of rigid disks. In simulations or experiments in 
2D geometry, a weak polydispersity is necessary to avoid 
local crystalline order. In our theoretical description, as 
we shall see below, the geometrical disorder is explicitly 
taken into account through the coordination number and 
the contact directions. We may thus assume that parti- 
cles are of the same size, and the steric hindrance will be 
characterized by a single exclusion angle 59 — tt/3. 

The directional organization of the contact network is 
described by the probability density of contact normals 
p{9) 0-Q]. This is the probability that a given particle 
has a contact along the direction parametrized by 6, the 
polar angle of the contact normal H; figure |l|. This dis- 
tribution is induced by the relative motions of particles, 
and so it essentially reflects the deformation history or 
the dynamics of the preparation process. The function 
p(9) is 7r-periodic. It is often assumed that it can be 
approximated by its truncated Fourier transform j(| 

p(9) = -{l + acos2(9~9 p )}, (1) 

where a represents the anisotropy of the texture and 9 P 
is the average contact direction. This functional form 
is reasonable for a simple deformation history, but it 
may considerably deviate from this simple form otherwise 
In order to bypass a particular fitting form such 
as equation |l| for characterizing the anisotropy, one may 
more generally construct the "fabric tensor" F = (n (g> n) 
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where the brackets denote averaging over all contacts in 
a representative element of volume, and (g) is the dyadic 
product Then, the texture anisotropy a is defined 

to be a — 2\Fi — F2I where F\ and F2 are the eigenvalues 
of F. 

A richer description of the texture is provided by the 
probability density f z (9) that a particle with z contacts 
has one contact along 9. This function is defined over the 
range [0, 2tt\. Like p(9), it can be Fourier expanded, but 
the corresponding anisotropy coefficients now will depend 
on z. At leading order in 9, we have 

f z (6) = ^-{l + a z cos2(9-6 z )}. (2) 

It is important to point out here that the steric exclu- 
sions do not allow for arbitrary global distributions f z . 
For f z defined by Equation 0, the anisotropy a z could 
take any values in the range [0, 1] (a z > 1 implies neg- 
ative values of f z for some directions, whereas negative 
values can be avoided by turning 9 Z into 9 Z + n/2). Nev- 
ertheless, the steric exclusions impose an upper bound 
a m ax on the anisotropy. Indeed, at most one contact 
can be found within an angular sector of 7r/3 around a 
particle, i.e. 

/ zf z (9)d9 < 1 (3) 

which together with the expression of f z given Eq. (j^) 
yields 

a max (z) (4) 

This means that a max decreases with z and becomes zero 
for z = 6. In order to evaluate f z for a granular sample, 
the subset of particles with z contacts are to be singled 
out. Our numerical simulations show that a z indeed de- 
creases with z ]To[ |. 

None of the above functions however directly provides 
an exhaustive representation of the local environment of 
a particle. A complete description requires the coordina- 
tion number z of the particle and the directions 9\ , . . . , 9 Z 
of the contact normals around it. Those angles should 
fulfill the constraint of steric exclusion around particles, 
so that an interval of width 7r/3 centered on the normal 
direction of any contact is excluded. In this way, the par- 
ticle environments should be described through multicon- 
tact probability distributions g z (6i, ■ ■ ■ , Z ). This distri- 
bution is 27r-periodic and its integration over all contact 
directions but one, 9k, gives back the global distribution 
f z {0k) defined above: 

f 9 z ({6i})d{0i}^ k = f z (9 k ), (5) 

where £ z {k) is the domain corresponding to integration 
from to 2tt over all the 9i except 9 k . 



Here, we would like to construct the local distributions 
g z {{Si}) from the global distribution f z {9). Of course, 
the solution is not unique since g z contains a much richer 
information than f z . The point, however, is to get a solu- 
tion that incorportes the required local information with 
no bias towards a particular solution. In the absence of 
local steric constraints, the most unbiased situation im- 
plies that 

z 

gM, ■■■,o z ) = l[f z (e i ). (6) 

In our case, this solution is wrong since the steric exclu- 
sions require 

5z ({^})=0 if \9 l -9 ] \<ir/3, (7) 

for i ^ j which is obviously not satisfied by the above 
solution. 

The solution that we propose is still to resort to a 
similar least biased estimate, provided the above steric 
constraints are taken into account. Operationnally, this 
translates into the maximization of an entropy functional 
S[g z ] under constraints. According to Shannon's entropy, 
we have 

%,]=-/ g z m})\og[g z ({9 l })}d{9 l }, (8) 

where T> z — ([0, 2tt]) z is the integration domain. This 
functional should be maximized over the set of functions 
that fulfill both the steric exclusions, Eq. (|t]), and the 
normalization conditions, Eq. (pj). We note that the en- 
tropy formalism, as a basic statistical tool, has been pre- 
viously applied to granular materials in order to char- 
acterize other aspects of a random packing such as the 
coordination number and the void ratio |11| . 

Equation (0) implies that the maximization of S has 
to be restricted to the region A z of T> z that is allowed 
by steric exclusions. Indeed, the maximization applies 
only to the lacking information, and here the value of g z 
is already known (g z — 0) for the angles belonging to 
the excluded region B z — T> z — A z . Hence, the above 
restriction allows to account for the steric constraints. 

On the other hand, the normalization conditions (Eq 
|5|.) can be imposed through Lagrange multipliers. This 
leads to the maximization of the functional 

T[g z ] = S[g z ] - V P HO^Cilg,] d9 % (9) 

over A z , where the Xi(9i) are the Lagrange multipliers 
and 

Ci[g z ] = f z (9i) - [ g z {{0i})d{9 k } k ^. (10) 

The solution is determined by setting the functional 
derivative of T[g z ] to zero. Let us remark that here the 
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only correlations among the angles are the steric exclu- 
sions themselves. This means that the angles Oi can be 
considered as independant variables when they are re- 
stricted to the allowed domain A z . As a consequence, 
the solution is given by a product of identical functions 
h of a single angle in the allowed domain. In this way, 
the general solution is written as follows 




(11) 



fc=i 



where G(0) is a 27r-periodic function such that G(9) = 
for |0| < 7r/3, and G — 1 otherwise. These prefactors take 
care of the steric exclusions. Note that without these ex- 
clusions, we have G = 1 for all angles and the solution (||) 
is recovered with h z (9) = f z {9). 

The function h z is determined by the normalization 
condition (^) which becomes now the following implicit 
equation in terms of h: 



I f[h z (9 k )d{9 j } j7ii = f z (0 i ), 
J A.® fe=1 



(12) 



where A z (i) — £ z (i) H A z . This equation can be used 
for a numerical computation of h z for a given f z . We 
used a simple iterative scheme over the functions u n {9) 
satisfying 



u n+ x{9x) = eu n + (1 - e) 



fz{9x) 



Jam u n(^) • ■■u n (9 z )d9 2 ■■■d9 z 



(13) 



where e is a relaxation parameter. 

The solution h z of equation (|I^) is simply a fixed 
point of the latter recursion. Our iterative process con- 
verges smoothly to the solution when the latter exists 
and, moreover, its convergence domain matches exactly 
the set of admissibles values of the anisotropy according 
to equation (||) . We show one example of solution in Fig- 
ure]^: The function h^{9) was computed according to the 
above procedure with fy(9) given by Eq. (0) for z = 3 
and 03 = 0.3. We see that h% has the same monotony as 
/3, but it deviates from a simple sinusoidal form. 

Of course, it is desirable to compare the multicontact 
distributions g z predicted by the above approach with 
direct data from experimental or numerical observations. 
However, for a sufficient statistical precision requires very 
large samples. In fact, given a granular sample, the sub- 
set of particles with z contacts are considered. Assuming 
that in a 2D system of disks, the number of particles is 
the same for each of the prevailing values 3, 4 and 5 of 
z, which is actually not the case since there are typically 
more particles with four contacts, then only nearly one 
third of the particles is available for the evaluation of g z . 
On the other hand, subdividing the interval [0, 2ir] into 
n angular sectors and requiring on average m events in 



each elementary box in [0,27r] z , one would need a sam- 
ple of m(n z ) contacts, or approximately (2m/z)n z par- 
ticles. For z = 4, n = 20 and m = 100, this amounts to 
8.10 6 particles, which is practically inaccessible with the 
present-time computation power. 

Nervertheless, we still may get a reliable comparison 
by focussing on the distributions p,i(A0) of the differ- 
ence A0 = 1 0i — 02 1 between contact directions extracted 
from pair distributions P2(0i,02)- The latter was calcu- 
lated from numerical samples and from the theoretical 
estimation of g z . To do so, we performed numerical sim- 
ulations of a system of 4000 particles with bi-periodic 
boundary conditions by means of the molecular dynamics 
method. We used a viscous-regularized Coulomb friction 
law and a linear spring-dashpot model for particle inter- 
actions ||. The coefficient of friction was 0.5 and the 
particle radii were uniformly distributed between R\ and 
i?2 with R-2 = 2R\. The system was subjected to simple 
shearing and thanks to bi-periodic boundary conditions 
we were able to keep shearing in the steady state for very 
large strains without any wall effect. Figure || shows the 
function fi(9) corresponding to the subset of particles 
with four contacts. The global steady-state anisotropy, 
a z , calculated from the fabric tensor was about 0.13. 

The distribution pd was calculated over the cumulated 
data from several well-separated snapshots of the steady 
state for the set of particles with coordination number 
z. For the theoretical evaluation of pd corresponding to 
these particles, we used the numerical distributions f z (9) 
for the calculation of g z following our approach, from 
which the distributions pi and pd were extracted. 

In figures || and ||, we have plotted theoretical and nu- 
merical pd as a function of A0 for z = 4 and z = 5, 
respectively. Although only the geometrical constraints 
are taken into account for the theoretical construction of 
g z (no force considerations), we see that the predicted 
distributions pd both for z = 4 and z = 5 fit fairly well 
the numerical data. The fit is nevertheless better for 
z = 5 than for z = 4. On the other hand, the fit for 
2 = 3 (not shown) is much less satisfactory. In fact, in 
going from z = 3 to z = 5, the requirement of force 
balance becomes less and less stringent, whereas steric 
exclusions dominate more and more the behaviour. Let 
us also remark that for z = 4 a small peak appears on the 
numerical curve at A0 = 27r/3, as can be seen in figure 
[|. This peak, which is absent from the theoretical curve, 
is probably the signature of a local order induced by the 
rather small polydispersity of our sample. 

These results suggest that, if the steric exclusions are 
correctly taken into account, a satisfactory evaluation of 
the statistics of local environments in terms of the func- 
tions g z can be obtained at least for z = 4 and z = 5. The 
steric exclusions control to a large extent the local distri- 
butions (figures f|, ||) and reduce the range of admissible 
anisotropies (equation ||) |l(J . 

A natural extension of this model, at the cost of addi- 
tional parameters, is to include the force balance as well 
as the steric exclusions. This would allow to tackle the 
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case of z — 3 more accurately than we did in this paper. 
The idea is to include both contact directions and con- 
tact forces in the model. The forces are coupled to a local 
stress tensor in the same way as the contact directions 
were coupled to the fabric tensor. Then, the entropy 
of local distributions is maximized regarding constraints 
arising from steric exclusions, but also the force balance, 
and constitutive laws such as Coulomb friction and Sig- 
norini's condition (positivity of normal forces). 

Apart from its importance in itself as a general model 
for the granular texture including steric exclusions, the 
characterization of the texture in terms of local environ- 
ments is a crucial step in order to relate the macroscopic 
behaviour to the particle scale. For example, the range of 
admissible stresses can be studied by considering repre- 
sentative configurations generated according to the mul- 
ticontact statistics of local environments []1(J . A similar 
procedure based on entropy formalism may be applied as 
well to other relevant microscopic configurations such as 
the cells composed of contiguous particles that carry the 
most local strains at the microscopic strains. 
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FIG. 1. Schematic representation of a particle environment with a coordination number 2 = 4. 6i is the angular position of 
particle i and 88 is the exclusion angle 
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FIG. 5. Distribution p(A9) of the difference between successive contacts of a particle with five contacts as predicted by our 
model and obtained from numerical simulations. 







